function [] = cross_validate()
% a script to perform Gaussian Process Classification, by Mark Norrish
% see Rasmussen & Williams for algorithms' derivations
% performs cross-validation on a data set

%disp('Loading data and initialising stuff...');

D = load('gp_data'); 
D = D(randperm(size(D,1)),:); % shuffle it to randomise trials
folds = 5; % or whatever
dim = size(D); 
size_fold = floor(dim(1)/folds); 
n_learn = (folds-1)*size_fold; % test data
ntot = n_learn+size_fold; % total data
D = D(1:ntot,:); % throw away data if need be
dim = dim(2)-1; % i.e. data are from a dim-dimensional space; the -1 is not counting the class label
n_class = max(D(:,dim+1));
func = @covSEard; % restricted for test purposes
Hyps = load('gp_hyps');
correct = 0;

approxF = zeros(n_learn*n_class, 1);

for fold_number = 1:folds
  X = [D(1:(fold_number-1)*size_fold,1:dim);D(1+fold_number*size_fold:ntot,1:dim)];
  
  Y1 = [D(1:(fold_number-1)*size_fold,1+dim);D(1+fold_number*size_fold:ntot,dim+1)];
  y = reshape(kron(ones(n_learn,1),1:n_class),n_learn*n_class,1)==repmat(Y1,n_class,1); % expressed as 1-of-n_class encoding
  % trust me, it works
  
  test_range = (fold_number-1)*size_fold+1:fold_number*size_fold;
  xs_star = D(test_range,1:dim); % test points
  ys_star = D(test_range,1+dim); % their labels
  
  K = zeros(n_learn,n_learn,n_class);
  ks_star = zeros(n_learn, n_class, size_fold);
  
  % this can't be further vectorised, can it?
  sigma_noise = 1e-5;
  for c = 1:n_class
    K(:,:,c) = func(Hyps(:,c), X, X) + sigma_noise*eye(n_learn);
    ks_star(:,c,:) = func(Hyps(:,c), X, xs_star);
  end  
  
  %fprintf('Computing f and pi of iteration #%i...\n', fold_number);
  [approxF, pi, E, M] = alg_3_3(n_learn, n_class, K, y);
  % time for Alg. 3.4
  % note lines 1-7 are already done; this is whence the test/test:= construct
  % in alg_3_3, rather than a generic for loop
  
  correct = correct + alg_3_4(E, func, Hyps, ks_star, M, n_learn, n_class, reshape(pi,n_learn,n_class), size_fold, xs_star, reshape(y,n_learn,n_class), ys_star);
end
%fprintf('I made %i correct predictions, out of a total of %i.\n', correct, ntot);
%fprintf('This is %f percent accuracy.\n', 100*correct/ntot);
dlmwrite('cv_results', [correct, ntot]);